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lUTKODUCmlUTY O^HE 

OHIGINAL PAGh A LANOSAT DIGITAL IHAGE RECTIFICATION STSTEM 


Peter Van Wie and Maurice Stein 


Goddard Space Flight Center 
Greenbelt, Maryland 


I. ABSTRACT 

DIRS Is a Digital Image Rectification System 
for the geometric correction of Landsat Multi- 
spectral Scanner digital Image data. DIRS removes 
spatial distortions from the data and brings It 
Into conformance with the Universal Transverse 
Mercator (UTM) map projection. Scene data In the 
form of landmarks or Ground Control Points (GCPs) 
are used to drive the geometric correction algor- 
ithms. The system offers extensive capabilities 
for “shade printing" lo aid In the determination 
of GCPs. Affine, two dimensional least squares 
polynomlnal and spacecraft attitude modeling tech- 
niques for geometric mapping are provided. Entire 
scenes or selected quadrilaterals may be rectified. 
Resampling through nearest neighbor or cubic con- 
volution at user designated Intervals Is available. 
The output products are in the form of digital tape 
In band Interleaved, single band or CCT format In 
a rotated UTM projection. The system was designed 
and Implemented on large scale IBM 360 computers 
with at least 300-500K bytes of memory for user 
application programs and five nine track tapes plus 
direct access storage. 


II. INTRODUCTION 

A. Overview 

The need for geometrically corrected Landsat 
MSS digital data Is strongly felt In a nuiriber of 
remote sensing application areas. A study of the 
1973 and 1974 Mississippi River floods motivated 
the development of the Digital Image Rectification 
System at Goddard Space Flight Center.* This 
system has gone through several stages of develop- 
ment and Is currently operational at Goddard. It 
is also offered for sale through COSMIC, the NASA 
software distribution facility at the University 
of Georgia. This paper will present the approach 
used In the development of DIRS, the distortion 
sources present In the data, the techniques which 
have been Implemented, and the results achieved 
thus far. 


B, Approach 

The approach used In the development of DIRS 
was pragmatic. Me were attempting to create a 
system which met our own needs within the context 
of available computing equipment. Our computer 
environment consists of two large scale IBM 360s: 
a model 91 and a model 75. Each has two mega- 
bytes of memory and a large collection of peri- 
pherals. 

Image rectification, as opposed to registra- 
tion, was desired since It allows Image data to be 
combined with map data. Combinations such as 
photographic overlays were originally contemplated, 
but direct digital combination of map and Image 
data was soon recognized as an Important possi- 
bility. Tne rectification capability also allowed 
for multi-scene registration and thus was selected 
as the preferred approach to geometric correction 
of the MSS Image data. 

The first choice to be made was that of a map 
projection. The Universal Transverse Mercator 
(UTM) map projection was selected as the most 
appropriate. The UTM projection offers the advan- 
tages of a cartesian reference system with metric 
measure and Is the primary projection for U.S.G.S. 
paper and digital topographic maps. See* for 
details on the UTM projection. 

with these considerations In mind, we col- 
lected an assortment of existing techniques and 
Implemented them. These techniques were derived 
from various sources Including work done for f<ASA 
by IBM’, TRU", CSC’, as well as other routine 
photogrametric techniques and work done by the 
Defense Mapping Agency’ •*•*. The development 
of new techniques was not our objective and was 
only employed where absolutely necessary. 

C. Objecti ves 

Our major objectives were: 

1. To rectify full or partial Landsat MSS scenes 
with the maximum attainable accuracy and produce 
digital output products suitable for further 
machine processing and analysis. 
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2. To use Available computing hardMare, support- 
ing software, and Image processing techniques. 

3. To produce efficient software which can pro- 
cess a full Landsat scene (30 megabytes of data) 

In an acceptable length of time. 

4 . To provide adequate flexibility to the user 
Including the ability to trade off GCP location 
effort for accuracy to allow them to achieve the 
level of accuracy they require at the lowest cost. 

III. DISTORTION SOURCES 

Distortions exist In the MSS digital Image 
data due to the combined effects of sensor opera- 
tion, orbit and attitude anomalies. Earth rotation 
and from atmospheric and terrain effects. No 
attempt Is made In DIRS to correct for atmospheric 
or terrain Induced geometric distortions.’**' 

A. Sensor Operation 

1. Aspect Ratio. An aspect ratio of approxl- 
mately 1.4:1 Is Introduced by the along scan 
oversampling of the MSS. The Instantaneous field 
of view (iFOV) Is 79 x 79 meters but along scan 
pixel separation (center to center) Is only 57 
meters. The line to line pixel separation Is 79 
meters giving no overlap or gapping In the cross 
scan direction. 

2. Mirror Velocity Nonlinearity. The MSS oscll- 
latlng scaning mirror used to sweep out the Image 
lines does not maintain a precisely constant 
angular velocity. Therefore, along scan distor- 
tions are introduced causing pixel compression 
or pixel stretching at various places long the 
scan line. This Is perhaps the most difficult 
distortion to correct since It Is poorly measured 
and changes gradually with time. Two standard 
mirror models are selectable In DIRS and provis- 
ions art! made for user defined mirror models. 

3. B^d to Band Misregistration . The MSS senses 
six scan~TTnes In each of four spectral bands on 
each mirror sweep. The light reflected off the 
scanning mirror falls on the end of a 4 x 6 array 
of fiber optic tubes which conduct it to the 
spectral filters and detectors. Displacement 
along the four element axis of the fiber optic 
array Induces an along scan offset In the observed 
ground position. This spectral misregistration Is 
corrected (to within a pixel) by the Introduction 
of pad pixels by the NASA Data Processing Facility 
(NDPF) where CCTs are generated from the video 
data tapes. A subpixel misregistration remains 
which Is not corrected by DIRS. 

4. Sensor Delay. The 24 MSS detectors are 
strobed sequentially by a precisely timed clock 
to Initiate the A/D conversion of the radiance 
data. The small time delays between the readout 
of each of the six sensors In a given band allow 
the scanning mirror to move forward a small 
amount. This gives rise to a small along scan 
displacement from the first through the sixth 
scan line In each swath. 


B. Orbit and Attitude Anomalies 

Significant distortions are caused by Incon- 
sistencies In spacecraft orbit and attitude. The 
major component of these distortions Is linear 
and can be co-rected with a single affine trans- 
formation of the entire Image using three ground 
control points (GCPs). However, small continuous 
variations In altitude, velocity, yaw, pitch and 
roll add nonlinear distortions on the order of 
20D to 300 meters over the scene. 

C. Earth Rotation 

Image Skew . Af the Landsat spacecraft 
travels southward It scans the Earth from West to 
East. The rotation ot the Earth causes each 
successive mirror sweep to begin a bit farther 
to the West. The overall geometric effect Is to 
skew the Image. The magnitude of this distortion 
Is proportional to the cosine of the latitude and 
Is greatest at the equator. 

Rotational Delay. The scanning of six lines 
with each mirror swMp causes the Earth rotational 
skew of an Image to be a step function. The Earth 
rotation correction In DIRS Is accomplished In 
two steps. The major portion Is corrected by 
treating the distortion as a linear (continuous) 
function over tbe entire scene area. This skews 
the Image as described previously. This linear 
correction leaves a residual saw tooth distortion 
which Is corrected (along with sensor delay 
distortion) when the Image resampling Is done. 

D. Miscellaneous 

Other geometric adjustments are necessary to 
compensate for the effects of Earth curvature, MSS 
versus map perspective and projection. 

E. Distortion Categories 

The distortions present In digital MSS data 
can be grouped according to the way In which they 
effect Image geometry. Three general categories 
can be defined. 

Global Continuous. This category Includes all 
those sources which operate over the full extent 
of a scene and are mathematically continuous. 
Attitude, orbit and aspect ratio are examples of 
distortion sources which fall in this category. 

Swath Continuous. These are distortions 
which operate In a consistent and continuous 
manner on each set of six scan lines In a mirror 
swath. Examples are mirror velocity nonlinearity 
ar.d Earth curvature. 

Swath Discontinuous . These are distortions 
which are discontinuous In nature and operate In 
a regular pattern within the lines of a swath or 
between two adjacent swaths. Sensor delay and 
the residual Earth rotation distortion (rotational 
delay) are In this category. 
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IV. TECHNIQUE 

The recti ficeliun proce&s used In DIRS 
consists of six najor steps end a number of special 
techniques. 

A. Location of Ground Control Points 

Ground control points are features or land* 
marks which are visible In the image and whose map 
coordinates can be determined. They are used to 
relate Image geometry to the desired map projec- 
tion. The Image coordinates of a landmark are the 
sample and line number (I.e., column and row) In 
the digital Image array where the landmark Is 
located. Hanual and automatic techniques for 
determining GCP Image coordinates are provided 
In DIRS. 

The manual method Involves the generation of 
shade prints of the GCP area with an Image coordin- 
ate border. These shade prints are made on the 
high speed line printer and use a combination of 
normal print chartcters to produce 16 levels of 
gray. Shade prints can be generated at full re- 
solution (one print position per pixel) or at 
expanded resolution (3x3. 5x5, 10x10 print posi- 
tions per pixel) using cubic convolution inter- 
polation. The edges of the snade print are 
annotated with the correct sample and line co- 
ordinates so that If a feature can be seen on 
the shade print Its coordinates can be measured. 

The automatic method of determining GCP Image 
coordinates Is accomplished through edge correla- 
tion. The GCP feature must, first be Identified 
using the manual technique and then extracted and 
entered on the GCP library tape. A second scene 
over the sar« ground location can then be processed 
using the automatic GCP location caoablllty. 

The edge correlation technique was developed 
for NASA by Computer Sciences Corporation for the 
Large Area Crop Inventory Experiment (LACIE) and 
Is a very '•ellable technique," It is based on the 
generally valid assumption that Image edges are 
Invariant with time. It Is well known that other 
Image qualities such as color. Intensity, or 
texture can change dramatically with time and 
season. These variations Introduce much of the 
difficulty In using other correlation techniques. 
Since edges represent the shape of ground features 
they are relatively Invariant with time. 

The edge extraction process begins with the 
computation of a reflectance gradient value at each 
pixel location. This Is done over the target sub- 
image area (from the manual GCP location) and the 
search area In the second scene. A gradient histo- 
gram Is then constructed for each of the subareas. 
The user specifies a threshold t determine the 
percentage of edge pixels. Genei-ally, about twenty 
percent of the pixels are used as edge pixels. The 
gradient histogram Is then analyzed to find the 
nearest breax point to split the gradient range In- 
to the desired percent of edge pixels and non-edge 
pixels. This break point will be different In the 
target area histogram than It Is in the search area 
histogram but will produce approximately the same 


percentage of edge pixels In both areas. The 
pixels whose gradient value exceeds the break point 
value are then coded as blnarv ones and all other 
pixels are binary zeros. The resulting blnarj edge 
Images are the. cross correlated and the point of 
maximum correlation Identified. The GCP coordin- 
ates In the target subimage are then mapped Into 
the search area and this point represents the image 
coordinates of the GCP In the second scene. Ac- 
curacies of plus and minus one sample and line are 
achieved with edge correlation. The reliability 
of the method Is In excess of eighty percent. 

The final operation In the location of GCPs 
is the determination of nap coordinates. General- 
ly seven and a half minute U.S.G.S. topographic 
maps are used fnr this purpose since they are very 
accurate and provide UTM coordinate tick marks. 

For maps without UTH coordinates the latitude and 
longitude values may be used. DIRS converts these 
geodetic coordinates to UTH automatically. Care 
must be taken to Insure that the same part of the 
GCP feature Is used for determining both map and 
Image coordinates. For example. If a road Inter- 
section Is used, the centroid of the Intersection 
Is usually selected as the precise GCP location 
and the coordinates of this centroid are determined 
In both the map and Image. 

B. Create Global Happing Functions 

The ground control point data developed In 
step A Is used to define functions which make 
transformations between the map and Image spaces. 
Three types of global mapping functions are avail- 
able In DIRS, affine transformation, two dimen- 
sional least squares polynomials and attitude 
model functions. The choice of the most appro- 
priate method depends on the ri.iiiber of GCPs and 
their distribution over the Image area. 

Affine Tra nsfor mation . The affine transfor- 
matlon is one wKTcTTmaps a triangle from one two 
dimensional space Into a triangle In a second 
two dimensional space. The mapping Is con- 
structed such that the verticles of the first 
triangle map Into the verticles of the second 
triangle and all other points within, on and 
outside the first triangle map Into propcrtlonata 
locations relative to the second triangle. An 
affine transformation accounts for the following 
distortions : 

Translation 
scale change 
rotation 
aspect ratio 
skew 

The mapping Is linear in the sense that straight 
lines map Into straight lines. Hathematical 1y . It 
Is not a true linear function since the translation 
term violates the linearity condition 

f(atb) • f(a) ♦ f(b) (I) 

Affine transformations are expressed as functions 
of the form 
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X • X„ ♦ i.u ♦ b,V 
0 1 1 


( 2 ) 

y - ♦ igU ♦ b2V (3) 

Mhert (X.Y) represent coordinates In one space and 
(U.V) coordinates In the other space. A set of 
affine triangles formed by the ground control 
points produces a pleceMise-llnear global mapping 
function over most of the Image area. Borders out* 
side the netMork of affine triangles are mapped 
using a large "master* triangle formed by three 
widely spaced GCPs. 

Some limitations of the affine transformation 
are that It Is a piecewise linear approximation to 
the "true* nonlinear distortion function. Also, 
the ground control points required to form an 
adequate set of affine triangles can be difficult 
to obtain. Some large Image areas may be sparce 
In features suitable for GCP usage. The affine 
method Is preferred where either small Image areas 
are needed from the rectified scene or where high 
accuracy Is not mandatory. 

Two dimensional Least Squares . In this method 
four functions are computed using a least squares 
fit to the GCP data. These functions are 


S • f^ (E.N) (4) 
L - f2 (E.N) (5) 
E ■ fj (S.L) (6) 
N - % (S.L) (7) 


where S Is the sample coordinate, L the lire 
coordinate, N the Northing coordinate and E the 
Easting coordinate. These functions are poly- 
nomials In two variables and may be defined as 
1st, 2nd. 3rd, 4th or 5th degree polynomials. The 
number of GCPs available determines the maximum 
degree of the polynomials. The polynomials have 
the general form: 

z - Cq ♦ c,x ♦ C2Y ♦ CjX^ ♦ c^xy ♦ CjY^ ♦ 

3 2 2 3 

CgX^ ♦ c^x^y ♦ CqXy‘ ♦ Cgy* ♦ . . . 

These functions are continuous and global and 
therefore offer advantages over the affine trans- 
formation If an entire scene or a large portion 
of a scene Is to be rectified. However, they have 
Inherent limitations due to the need for a large 
number of ground control points and tne require- 
ment that the edges and corners of an Image be 
well covered by GCPs. Some of the edge problems 
can be minimized by using "auxiliary" control 
points. ACPs are pseudo ground control points 
created at Image locations where GCPs are needed. 
The Image coordinates for an ACP and a suitable 
affine triangle (formed by existing GCPs) are 
specified by the user and the UTM coordinates 
at the ACP a'‘e computed using the affine trans- 
formation. ACPs are somewhat Inaccurate but with 
judicious placement they can constrain the least 
squares functions and Introduce very little error. 
Typical locations for ACPs are directly out from 


the Image corners at a distance of 500 sait^les 
and 500 lines. Some experimentation may be needed 
In a particular situation to achieve the best 
results. The remaining problem with the least 
squares approach Is the number of GCPs required 
to obtain the best results. In excess of thirty 
GCPs are needed to get good results with fifth 
degree polynomials. This problem is reduced 
through the use of automatic GCP location using 
edge correlation where possible. 

Attitude Model . The final option for a globa 
mapping function Is the use of an attitude model 
for the Landsat spacecraft. This Involves the 
determination of functions relating time to space- 
craft position, velocity, altitude, yaw, pitch and 
roll. DIRS has the capability to use these func- 
tions to geometrically correct the Image data, but 
It does not have the ability to compute the func- 
tions from the GCP data,*' The Goddard Trajectory 
Determination System (GTDS)*' Is used to compute 
these functions from the GCP data. Approximately 
nineteen GCPs are needed to obtain the best fit 
for all of the attitude and orbit functions. This 
method offers the advantages of requiring fewer 
GCPs and having a somewhat lower sensitivity to tht 
GCP dispersion pattern. 

C. Create an Interpolation Grid 

The global mapping functions computed In step 
B could be used directly to transform the 1it<age 
data Into the desired map projection. This Is not 
practical, however, due to the enormous amount of 
computation required to evaluate these functions 
over millions of points. The solution to this 
problem Is the creation of an Interpolation grid. 
The functions can be evaluated at each mesh point 
In a grid covering the Image area and then an 
efficient bilinear Interpolation technique can 
be used to map points within each grid cell. 

The Interpolation grid Is composed of approxi- 
mately twenty horizontal and twenty vertical grid 
lines. These lines are defined In the UTH map 
space by functions of the form 

y ■ MX ♦ B (7) 

and AX « BY » C • Q (8) 

where Y Is Northing and X Is Easting. Tne UTM 
Intersection coordinates for each pair of norlzon- 
tal and vertical grid lines Is then computed. These 
Intersections (cr mesh points) are then transformed 
through global mapping functions to compute the 
corresponding Image coordinates. When each mesh 
point Is known In both UTM and Image coordinates 
we are ready to map the Image data in an efficient 
manner. The mapping Is actually an Inverse map- 
ping, I.e. It proceeds from if, e map space (output 
coordinate system) back Into the Image space (In- 
put coordinate system). This arrangement has the 
advantage of allowing Interpolation to occur In 
the Input (Image) space Instead of the output (map) 
space. This proves to be a much more convenient 
approach than forward mapping. 
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0. Step through Interpolation Grid 

The desired end product of < digital Image 
rectification Is ainlformly spaced matrix of 
pixels which produces an Image In confoneance with 
a (UTN) map projection. It Is Impossible to main- 
tain a uniform lattice in both the Image and map 
space since the mapping between the two Is non- 
linear. The regular lattice of points In the map 
space form a nonuniform lattice In the Image 
space. Trie objective of stepping through the 
Interpolation grid Is to identify the map co- 
ordinates of the uniformly spaced points In the 
map space. These are the locations at which Image 
Intensity values are desired for the rectified 
output Image. 

E. Compute Image Coordinates 

The map location selected In step D Is now 
transformed Into Image coordinates using bilinear 
Interpolation within the local grid cell. Pad 
pixels are Introduced at the beginning and end of 
each resample line to form a rectangular output 
pixel array. The Image coordinate computation Is 
done by breaking the resample line Into segments. 
One segment corresponds to the portion of the 
resample line between two vertical grid lines. 

The Image coordinates at the Intersection of the 
resample line and the two vertical grid lines are 
computed and an equation of the line segment In 
the Image space Is derived. Actually, two equa- 
tions are derived, one for the sample coordinate 
and one for the line coordinate. These equations 
are evaluated at each resample location to produce 
the Image coordinates at which resampling Is to be 
done. 

F. Resample 

Resampling refers to the determination of an 
Image Intensity value at a given location. 
Typically, this location falls between and not on 
exact pixel centers and some form of Interpolation 
Is required. Two Interpolation techniques are 
available In DIRS. The simplest and most efficient 
method Is Nearest Neighbor. In this method the 
pixel whose center Is nearest to the resample loca- 
tion Is used to supply the Intensity value at the 
resample location. This method Introduces up to 
one half sample and line of geometric error. For 
a large subimage area or full scene this Is 
acceptable and produces adequate results In a 
short period of computation time. 

The second resampling technique provided In 
DIRS Is cubic convolution. This method was de- 
veloped by TRW and Is an efficient approximation 
to the theoretically optimum Interpolation using 
s1n(x)/x. While cubic convolution Is much more 
efficient than s1x{x)/x It Is a great deal slower 
than nearest neighbor. 

Whether resampling Is done using nearest 
neighbor or cubic convolution, DIRS allows the 
option of setting the spacing between pixels and 
lines. Since Landsat pixels are spaced 57 x 79 
meters and since even spacing Is desired In the 


output array. It Is necessary to Introduce some 
redundancy to preserve all the Information In the 
original Image. A sample and line spacing of SO 
meters Is frequently selected for output pixels. 
Fifth meter spacing gives a scale of 1;1, 000,000 
when the Image data Is displayed on a film re- 
corder with a SO micron toot size. Soaclnq at 
small as IC .Tcters has been used successfully with 
DIRS for small subimage areas. 

6. Special Techniques 

Steps A through F above describe the general 
approach used In DIRS. This covers the correction 
of all distortions In the category ‘global con- 
tlnuous'as discussed In Section III-E. Special 
methods are Implemented to remove swath continuous 
and swath discontinuous distortions. 

Since swath continuous error'> could also be 
viewed as global continuous errors they could 
be corrected through the global mapping functions. 
However, because these distortions are highly con- 
sistent from swath to swath, there Is no need to 
"load" the global mapping functions to correct 
them. Whenever distortions can be removed without 
resorting to scene data (GCPs) It Is advantageous 
to do so. This reduces the total number of GCPs 
needed to achieve the same accuracy. The method 
used for swath continuous distortion corrections Is 

1. Adjust the sample coordinate of each 
GCP by subtracting the swath continuous 
errors. 

2. Generate the Interpolation grid usin.., 
the adjusted GCP data. The Image co- 
ordinates at mesh points In this grid 
now correspond to an Image which does 
not have swath continuous distortions. 

3. Reintroduce the swath distortion by 
adding It to the sample coordinate at 
each mesh point In the Interpolation 
grid. This forces the grid to map 
back Into the correct (raw) Image 
coordinates In the real Image space. 

The effect of these operations Is to super- 
impose the swath continuous distortion functions 
on the global mapping functions. 

The thirr’ category described In Section III-E 
Is swath discontinuous distortions. These distor- 
tions must be corrected In the resampling process. 
Both nearest neighbor and cubic convolution re- 
sampling algorithms In DIRS correct frr sensor 
delay and Earth rotational delay distortions. The 
technique used Involves reposi tlonlr j the Input 
Image pixels prior to resampling. Thfl reposition- 
ing Is the Inverse of the original distortion and 
removes the distortion, but It also creates an Ir- 
regular array for the resampling algorithm. See 
Figures 1 and 2 to compare the original and re- 
post ;.1:ned pixel arrays used for cubic convolution 
resampling. 



An tddltloncl tcchn1()u« provldvd In DIRS U 
th« rcfflovdl of ttnsor deity tnd Etrth rotatlontl 
delay distortions from the expanded shade prints. 
This produces more realistic appatrlmj features on 
the shade prints and helps In determining the 
Image coordinates of the GCP features. This pro- 
cedure Introduces a small error In the Image co- 
ordinate associated elth the given feature. This 
error Is eliminated when the GCP table Is pre- 
processed prior to developing the global mapping 
functions. The sample coordinates of the GCPs 
are adjusted to put back the error due to sensor 
delay and Earth rotational delay which had been 
taken out of the shade prints. 

V. RESULTS 

DIRS has been used to rectify eight scenes 
thus far. Both Landsat 1 and 2 scenes have been 
processed. The accuracy of the results depend 
primarily on three factors: 

(1) The number of GCPs 

(2| The distribution of GCPs 

(3) The mirror model used. 
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VI. CONCLUSIONS 

DIRS Is operational and performing up to 
expectations. The objectives described In Section 
II-C have been achieved. Further work will be 
required to determine the limits of accuracy pos- 
sible with DIRS. Additional work could also be 
done to Improve efficiency but the existing run 
times are within acceptable limits. A full Land- 
sat scene can be rectified using nearest neighbor 
resampling at 50 meter steps In less than 5 min- 
utes of CPU and 1.6 minutes of I/O time on a 
360/91 comouter. This does not Include the run 
time for reformatting the four CCTs and generating 
shade prints. 
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